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In this paper we calculate the interfacial resistances to heat and mass transfer through a liquid- 
vapor interface in a binary mixture. We use two methods, the direct calculation from the actual non- 
equilibrium solution and integral relations, derived earlier. We verify, that integral relations, being a 
relatively faster and cheaper method, indeed gives the same results as the direct processing of a non- 
equilibrium solution. Furthermore we compare the absolute values of the interfacial resistances with 
the ones obtained from kinetic theory. Matching the diagonal resistances for the binary mixture we 
find that kinetic theory underestimates the cross coefficients. The heat of transfer is as a consequence 
correspondingly larger. 
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A number of different methods have been used to obtain the surface transfer coefficients for one-component systems: 
experiments plM^ . molecular dynamic simulations [25l - l29| , kinetic theory [30l - l33j . In a paper coauthored by one 
of us [2H the interfacial transfer coefficients were calculated with the square gradient theory for a one-component 
system, and compared to the data in the above references. Even for one-component systems the database of interfacial 
transfer coefficients is poor and these data are pretty scattered. The situation is even worse for mixtures. There are 
only few experiments available [H, [24{ at a very restrictive range of conditions, i.e. for instance, at infinite dilution. 
^ \ No molecular dynamic simulations are available yet. The only source of values of interfacial coefficients is kinetic 
. theory [32|, |33[. This theory is most appropriate for short range potentials and low density gases. There is evidence 
from molecular dynamic simulations for one-component systems for longer range potentials [29j that the coupling 
transfer resistivities for liquid-vapor interfaces of real fluids are substantially larger than those predicted by kinetic 
theory. 

It is the aim of this article to determine the heat and mass transfer resistances of the interfacial region. The values 
^ ' of these transfer coefficients, or even their order of magnitude, are extremely important for industrial processes which 
involve evaporation and/or condensation of mixtures. Among these processes is, for instance, distillation, when one 
needs to separate components with different volatilities. As this involves evaporation and/or condensation repeatedly 
many times, it is very important to know the exact effect of the surface. Some values of the interfacial transfer 
' . . coefficients may favor transport of a component, while other values may not. Of particular interest are the values of 
the cross coefficients, which contribute to reversible transport, and which are in most descriptions neglected jUj]. 



We will verify that integral relations, derived in [36| give the same values of resistances, obtained directly from a 
non-equilibrium numerical solution. The numerical solution is obtained using the non-equilibrium square gradient 



model 37[. It is desirable to compare our predictions with other methods, in particular molecular simulations and 
experiments. Such data are not available yet, however, and we will therefore use the predictions of kinetic theory to 
Jv>( \ compare with. 

In our approach we use the local resistivity profiles. The values of the local resistivities in the liquid and the 
vapor phases are chosen on the basis of experimental values. In the interfacial region there are small peaks in these 
resistivities. The results of molecular dynamics simulations (28j support the existence of such peaks in the local 
resistivities in the interfacial region. The amplitudes, being the adjustable parameters, control the magnitude of these 
peaks. The square gradient approach gives a natural tool to incorporate these peak in the theory. Possible values 
of these amplitudes are found by matching the diagonal transfer coefficients to values predicted by kinetic theory 
Using these amplitudes we find that the value of the cross resistivities is 1-2 orders of magnitude higher then the one 
from kinetic theory. The results indicate that kinetic theory underestimates the interfacial transfer coefficients in real 
fluids. One of them even has a different sign. 

Consider a planar interface between a liquid and a vapor of a mixture through which there is evaporation or 
condensation. The mixture is in a box with gravity g directed along the x-axis from left to right. The gas phase is 
therefore in the left part of the box and the liquid is in the right part. Due to evaporation or condensation there exists 
a mass flux of component i, which is equal to the mass of component i transferred through a unit surface area per 
unit of time. Furthermore, there exists the total energy flux J e , which is defined similarly. In stationary conditions 
these fluxes are constants (independent of x). 
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A surface can be described by Gibbs excess properties. We refer to [36| and [37} for an explanation how these 
quantities can be introduced in non-equilibrium. Due to non-equilibrium conditions the temperature and the chemical 
potentials of the components are not the same in the liquid and in the gas phases. Let T l and T 9 be the extrapolated 
temperatures of these phases at the surface. The exact position of the dividing surface is irrelevant for the following 
analysis. Furthermore, let [i\ and /i 9 similarly be the extrapolated chemical potentials of the i-th component at the 
surface. 

The paper is organized as follows. In Sec. [II] we discuss the different forms of the excess entropy production of 
the interface and introduce interfacial resistances. Sec. gives the overview of the expressions for these coefficients 
predicted from kinetic theory. We further build a procedure to determine the actual values of these resistances directly 
from a non-equilibrium numerical solution in Sec. |IV| and from integral relations, which use only equilibrium profiles 
in Sec. [V]. We compare the predictions of all three methods in Sec. [VI] and discuss the results in Sec. [VII] . 



II. EXCESS ENTROPY PRODUCTION 



In a previous paper [3q | we have obtained the following relation for the excess entropy production for the Gibbs 
surface in case of transport through the interfac^ll 



e \ T l Ta I ^ ^ \ Tl Ta 



rpi Ta 



(III) 



where fii = fii +v 2 /2 — gx s , with v the barycentric velocity and x s the position of the dividing surface. We introduce 
the measurable heat flux J' q by 



(11.2) 



i=l 



where hi = hi + v 2 /2 — gx s = Jli + Ts^, with Si the partial entropy and hi the partial enthalpy. Using the measurable 
heat flux on the vapor side, the excess entropy production can then be written as (36| 
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An alternative form of this expression is 
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It is important to realize that Eq. (|II.4|) . Eq. (|II.3|) and Eq. are exactly equivalent. It is common to do these trans- 

formations neglecting third and higher order contributions in the deviation from equilibrium. Such approximations 
were not needed here. If one neglects such higher order terms one may write Eq. (JO)) in the form 



U * ~ ( rpt Jig j ^2 Tl 
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(II.5) 



This expression is convenient if one wants to write the chemical forces in terms of the natural logarithm of the partial 
pressure divided by the partial vapor pressure of the liquid^. We refer to 351 for a discussion of this. 

Eq. (|II.5|) has the form of the entropy production for the surface used in [35J . It was obtained there using the local 
equilibrium hypothesis, which we have proven to be valid in [37| . In [36j | we have derived Eq. (]II.4[) independently, by 
calculating the excess of the continuous entropy production. 

We now consider a binary mixture. The excess entropy production can be written as 



u s = J^ 3 X q + J h Xf + J i2 X. 



(II, 



1 In |36H we have used the notation £ [cr s ] for the excess entropy production to distinguish it from the local entropy production <r s . Here 
we do not use the local entropy production and therefore will denote the excess entropy production by a s to simplify the notation. 

2 These partial pressures are defined as the molar concentrations times the total pressure. 
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where 

X n = 



1 1 



The resulting linear force-flux relations are 

X q = Rqq Jq 9 + Rql ^1 + Rq2 
'9 _ n5 I'.JXPS 7. _L R 9 



or in the matrix notation 









R 9 = ( 








we have 



(II.7) 



X( =R 9 lq J^ + B 9 nJ 5l +R 9 12 J S , (II. 



X 9 ee | Af | , R 9 ee ( i?| 41 i?f 2 | , J 9 ee ( J 6 I (IL9) 

R-2q R-21 "2 



X 9 =R 9 -J 9 (11.10) 

The resistance matrix R 9 satisfies the Onsager reciprocal relations, i.e. R ql = R 9 q , R q2 — R%^ an d R21 = ^12- 

In the above expressions for the entropy productions we used the measurable heat flux Eq. (III. 2 j) on the vapor side 
of the surface J q 9 ■ One can similarly use the measurable heat flux on the liquid side of the surface J q l - The resulting 
resistance matrix R £ differs from R 9 . We refer to [35[ for the details of the alternative procedure. 



III. KINETIC THEORY 

According to [HI, p. 180] kinetic theory gives the following expressions for the surface transport coefficients for a 
two component mixture 
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Wi = Aj/(Ai + A 2 ) 



where i? is the universal gas constant, Aj and cf are the thermal conductivity and the gas coexistence concentration 
of the i-th component respectively. <7j is the condensation coefficient of the z-th component, which arc parameters in 
this theory, and <5y is the Kroneker symbol. Furthermore, Mi, the molar mass of component i, appears in Eq. (IIII.2[) 
to adopt the molar transfer coefficients used in (35[ to the mass transfer coefficients used in this paper. All these 
quantities and as a consequence the resistances are calculated for a liquid and a vapor in coexistence at the temperature 
T and the chemical potential difference /Ui2 = /ii — ^2 of the surface, see (35l . p. 180] in this context. 
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IV. NON-EQUILIBRIUM CONTINUOUS SOLUTION 

Assume we have the numerical solution for a particular non-equilibrium stationary state. That is we know all 
the fluxes J 9 and forces X 9 used in Eq. (III.10[) : the constant fluxes are obtained directly from the non-equilibrium 
solution and the extrapolated bulk profiles are obtained using the procedure described in (37j . We now consider the 
following problem: to determine the transport coefficients for the whole surface having the non-equilibrium solution. 
This problem is, in a way, inverse to the common one, where one knows the resistances and, say, forces, and needs 
to determine the fluxes. As one can see, Eq. (|II.8|) has 9 unknown resistance^! and only 3 equations. It is therefore 
not possible to determine all the transport coefficients uniquely having only one stationary state solution. In order to 
incorporate more equations we need to consider other non-equilibrium stationary solutions which are independent of 
the previous. An important observation should be made here. 

In (3tJ we have verified the validity of the hypothesis of local equilibrium of the surface. This implies, that the 
resistance matrix R 9 is a function of thermodynamic parameters, say the temperature T and the chemical potential 
difference /Ui 2 , of the surface: R 9 = R 9 (T S , fJ,f 2 ). in [23I we saw, that the temperature of the surface and the chemical 
potential difference of the surface depend on both, the equilibrium temperature and the chemical potential difference, 
as well as on the size of the perturbation. Let 6 indicate the size of a non-equilibrium perturbatiorjj so that 

T" = T s (T eq , fi 12 ,eq;ft) (IV 1) 

Ml2 = Ml2( r eg,Ml2,eg;fi) 

Furthermore, X 9 = X 9 (6) and J 9 = J 9 (fi). In order to be able to use several independent perturbations as a source for 
the resistance coefficients, we must ensure that for all perturbations the temperature of the surface and the chemical 
potential of the surface are the same. The simplest way to ensure this is to assume that T s w T eq and fif 2 ~ /U12, eq- 
As is clear from Eq. (jIV.II) . this can be considered true if the perturbation rate fi is small enough. As we decrease 6, 
the accuracy of this assumption increases and in the limit B — >■ it becomes exact. It follows that 

R 9 = R 9 (T eq ,^ eq ) = limR 9 (T e(? ,/i 12 , e(? ;6) (IV.2) 

In practice there exists a particular size fi ei? of a perturbation, such that for all 6 < 6 e9 , T s w T eq and /if 2 ~ Mi2,eq 
with a satisfactory accuracy. 

One should also note that the accuracy of a particular numerical procedure may impose a lower bound for the size 
of the the perturbation 13 as well. All the non-equilibrium profiles and therefore forces and fluxes are calculated by 
solving the system of differential equations numerically with some particular accuracy. If a perturbation rate fi is 
lower then this accuracy, say B nMm , then the data obtained from the numerical procedure are not reliable. We may 
therefore use Eq. (|II.I0[) only if the perturbation rate 6 is in the range Q num < fi < fi e g- The boundaries of this range 
should be established empirically. 

We determine the transport coefficients from two different methods: a "perturbation cell" method^] and an 
experimental-like procedure. For ease of notation we will suppress the superscript g in the rest of this section, as the 
procedure is the same for vapor and liquid interfacial resistances. 

A. Perturbation cell 

Consider a stationary state which is perturbed from equilibrium by setting the temperature of the liquid^ T{x t ) — 
(I + f3r)T eq , the pressure of the gas p(x 9 ) = (I + (3 p )p eq and the mole fraction of the liquid ( e (x e ) = (1 + Pc)Ceq 
independently. The resulting non-equilibrium state is therefore a function of the parameters f3: 

X(Pt,P p ,P<) = R(T eg ,/i 12)eg )-J(/3 T ,/3 p ,/3 c ) (IV.3) 



3 Solving the inverse problem we have to ensure the validity of the Onsager reciprocal relations. This is one of the criteria to limit the 
size of the perturbation. This means that we have 9 independent resistances, but not 6. 

4 Note, that a non-equilibrium state can be achieved by perturbing several independent quantities simultaneously. In this case we have 
several perturbation parameters . . . , /3 P . A measure fi is a norm of this p-dimensional vector of perturbations. The exact expression 
for this norm is irrelevant, as soon as it goes to zero if and only if all ^1 , ... , 0p go to zero. 

5 This method was first used by Johannessen et. al. in I'i4l for one-component system. Here we discuss the grounds for the legitimacy of 
this procedure and generalize it to mixtures. 

6 One should not confuse T(x e ) with T e . The former is the actual temperature at x = x e , i.e. at the box boundary on the liquid side. 
The latter is the temperature extrapolated from the liquid phase to the interfacial region and calculated at x = x a , i.e. at the dividing 
surface. 
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where X, J and R are given by Eq. pi.9[) . Consider the following set of 8 independent non-equilibrium perturbations: 
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(IV.4) 



Consider now the 3x8 matrices X and 3 which contain 8 column vectors X and J respectively for each non-equilibrium 
perturbation specified above. For these perturbations X = X(B) and 3 = 3(P) are functions only of one parameter B. 
It follows from Eq. pVTI) that 

3Z(p)=R(T eq ,(, 12teq )-3(l3) (IV.5) 

where 8 should be in the appropriate range, as discussed above. In Appendix [X] we discuss the method to obtain 
this range. From Eq. PV.5P we obtain 

R(Teq,fH2,eq) = (X(B)^ T (l3))-(m^ T (P)y 1 (IV.6) 

where superscript T means the matrix transpose and _1 means the inverted matrix. 

We note, that in order to obtain the resistance matrix R uniquely, it is sufficient in principle to impose any 3 non- 
equilibrium perturbations which have sufficiently small perturbation parameters Bt, Bp and Be- This would give us 
3x3 = 9 independent equations. The method presented above makes the resistance matrix converge to R(T e9 , /ii2, e q) 
as fast as B 2 goes to zero, however. This is achieved by using 8 symmetric perturbations at the "corners" of a 
three-dimensional "perturbation cell", so changing B to — B does not change the "perturbation cell" and the resulting 
R. 

Because of using 8 perturbations instead of 3, there are 5 superfluous perturbations which make the system of 
equations (|IV.5I) to be overdetermined. Contracting both sides of Eq. (|IV.5|) with 3 T we actually average all the 
perturbations which are spread around T eq and fJ-i 2 ,eq in the least square sense. As the components of Z matrix are 
linearly independent, this guaranteers the matrix $-Z T to be invertible. Thus, the inverse matrix (Z'Z T ) _1 exists 
and Eq. (IIV.6|) is mathematically legitimate. In the numerical procedure the expression on the right hand side of 
Eq. pV.6p is obtained using Matlab matrix division /. 



B. Experiment-like procedure 

In experiments it is convenient to measure the corresponding coefficients by keeping zero mass fluxes through the 
system. It is also convenient to work with the total mass flux J m = + J^ 2 and the flux of one of the components 

= , rather then with fluxes of each component separateljQ, and J^ 2 . The excess entropy production Eq. pi.6[) 
can be therefore written as 

a s = J' q X q + J 5 A^ + J m X m (IV.7) 
where = X\ — Xi and X m = X2 . The resulting force- flux relations pi.lO[) have the following terms 

Rqq Rq£ Rqm 

X=\X S \, R = ( R iq R Sm I , J = I J e I (IV. 

, Rm.q Rm£ Rmm / 

where the resistances for different force definitions are related as 

Rq£~\~Rqm Rqm 
mq Rmm ~t~ Rm£ R^m Rmrn R^m J (IV.9) 
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7 One of the reasons for this is that it is hard to make only J* = 0, keeping Jj 2 finite. 
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Consider a stationary state which is perturbed from equilibrium by setting the temperature of the liquid T(x l ) — 
(f + P)T eq . The perturbation parameter /3 is a small number. The second perturbation constraint we impose is either 
= or C*( x ) = Ce ? an d we introduce the perturbation parameter which is in the former case and 1 in the 
latter one, which will be used as a subscript. The third perturbation condition is either J m = or p(x 9 ) — p eq and 
the corresponding perturbation parameter v m is or 1 respectively, which will be used as a subscript. The resulting 
non-equilibrium state is therefore a function of 3 parameters: 

X* tl „ m G0) = R(T eq , fJ,\2,eq)-^u 4 ,u m (l3) (IV. 10) 

where X, J and R are given by Eq. pV.8[) . 

Consider the following set of 3 independent non-equilibrium perturbations: 

Xqo(/3) = R(T eg ,^i2,e 9 )-Joo(/3) 

X10G8) =R(T e „/ 1 i2 >efl )-Jio(/3) (IV.ll) 

X n (/3) =R(Te 9 ,/ii2 >e9 )-Jii(/9) 

Further on for simplicity we will suppress arguments j3 and (T eq , \iyi, eq)- 
From the first of Eq. (|IV.11|) we find 

R qq = X q , 00 / Jg t 00 

R (q = X Ct oo /Jloo (IV. 12a) 

Rmq = V m 00 / Jq y 00 

From the second of Eq. pV.ll[) we find 
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(IV. 12b) 



The values X 10 and Jio are found directly from the calculations and the values of R qq , R^ q and R mq are those which 
are found in Eq. (|IV.12al) , given that the perturbation rate (3 is small enough. From the third of Eq. pV.ll[) we find 

Rqm = (-Xg, ll — R qq J q< n — R q £ n) / J m ^ n 

R$m — 11 — R$q Jq. \\ — R& J$, ll) / Jm, 11 (IV.12c) 

Rmm — \X m , 11 — Rmq Jq, 11 — Rm£ J$, ll) / Jm, 11 

Again, all the quantities on the right hand side of Eq. (IIV.12cp are known and we therefore can find the remaining 
resistivities. 



V. INTEGRAL RELATIONS 

In [HI we have established the general approach to derive integral relations between the surface resistances and 
local resistivity profiles. In this section we apply it to find the relations between the resistances R used in Eq. (|II.10|) 
and Eq. (|II.9[) and local resistivities r. Using the method described in [36[ we find 
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(VI) 
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where the operator l£ is defined as 



x 



£{(j>}{x s )= / dx[4>{x)- 4> 9 {x)Q{x s -x)-(j) e {s)e{x-x s )] (V.2) 



where <f> 9 and (j/ are extrapolated from the gas and liquid respectively profiles of <fi, while x 9 ' s and x l,s are the surface 
boundaries. 

This method requires the equilibrium profiles for the enthalpy h(x) and the mass fraction £(x) across the interface. 
Both of them could be easily obtained from the equilibrium square gradient model, see [37] for details. In contrast 
to the methods in Sec. |IV| . this requires calculating only the equilibrium profiles, but not the non-equilibrium ones, 
which is a much easier calculation. 

The integral relations also require the local resistivity profiles r qq (x), r q i(x), and ru(x) across the interface, which 
were modeled in the square gradient theory as 



(V.3) 



where qo(x) and q\{x) are modulatory curves for resistivity profiles which depend only on density profiles and their 
first derivatives. We refer for the details to [37]. qo(x) is a smooth arctan-like function which changes its value 
from to 1 within the range [x 9,s ; x e,s ] and qi(x) is zero on the boundaries of the [x 9 ' s ; x 1 ^] interval and has a peak 
proportional to the square gradient of the density inside this interval. Thus, the first two terms in each expression 
for the resistivity represents a smooth transitions from the gas bulk resistivity to the liquid bulk resistivity, while the 
third term represents a peak in the resistivity proportional to the square gradient of the density. For each resistivity 
profile r 9 and r l are the equilibrium coexistence resistivities of the gas and liquid phase respectively. They are related 
to the measurable transport coefficients such as heat conductivity, the diffusion coefficient and the Soret coefficient. 

The square gradient model used 3 adjustable parameters a qq , a q \, an which control the size of the peak in the 
resistivity profiles in the interfacial region. The interfacial resistance coefficients R will therefore depend on these 
coefficients, R = R(a gg ,a s i,aii), which we will investigate. 
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VI. RESULTS 

We consider a binary mixture of cyclohexane and n-hexane, as we did in [37j ■ 

We find in Appendix [X] that (3 = 2-10~ 4 is an optimum perturbation rate both in the "perturbation cell" and 
"experimental-like" methods, see Sec. [IV] . We have verified that both these methods lead to essentially the same 
values of the resistance coefficients. The numbers given below are taken from the "perturbation cell" method. 

Furthermore in Appendix [B] . we find the range of adjustable amplitudes ot qq , a\ q and an, for which the description 
is thermodynamically consistent. We find that a qq ~ 10, an ~ 1. The value of a.\ q is found to be irrelevant. 

In this section we suppress the superscript g for the resistances for ease of notation. 



A. Comparison to kinetic theory 

In this subsection we investigate the values of parameters a qq , ot\ qi an which makes the coefficients agree with 
the kinetic theory coefficients. We do it for j3 = 2e-4 as this perturbation rate gives the most accurate results. 
Furthermore we use the temperature T eq — 330 K and chemical potential difference /Ui2, eg — 700 J/mol. The values 
of parameters, used for kinetic theory are the same, as we use in our calculations. Particularly, the heat conductivities 
are Ai = 0.0140 W/(m K) and A 2 = 0.0157 W/(m K), M 1 = 84.162 g/mol and M 2 = 86.178 g/mol. 

We found that a variation of a\ q from to 10 makes the diagonal coefficients vary about 1 % and the cross 
coefficients vary not more then 5 %. As the variation of ai q is quite substantial, the variation in the coefficients which 
it induces is negligible. We therefore take a\ q = in all further analysis. 

Let us use subscript pc for the resistivity matrix obtained from the "perturbation cell" method and subscript kin 
for the resistivity matrix obtained from kinetic theory. For the above parameters R qq , pc = 2.96792 x 10 -11 . We found 
that Rqq^kin is practically independent on an while it depends linearly on a qq , see Fig. pQ. One can see from the 
plot, that they are the same for a qq ~ 9. 
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FIG. 2. Dependence of Ru, pc and Rii, V c on em (dots, bottom axes) and Ru,kin and R22,kin on cri and 02 (curve, top axes), 
respectively. Data are obtained at T eq = 330 and fj,\2, eq = 700 for a qq — 9 and ot\ q — 0. 



The diagonal coefficients i?n iPC and i?22,pc depend both on a qq and an. Since we have found the value of a qq 
already, we will further investigate the dependence of Rn, pc and i?22,pe using this value of a qq and varying only 
an. The diagonal coefficients i?n,fcm and R%2,Hn depend, in their turn, on the condensation coefficients u\ and 02 
respectively. We plot this dependence in the same plot with the dependency of Ra,p c (i = 1; 2) on an, see Fig. [2], 
The dependence of Ru, pc on an is given by the dotted line with the values of an drawn on the bottom x-axes. The 
dependence of Ru.kin on <Ji is given by the solid line with the values <ii drawn on the top x-axes. 

Consider a particular value Ru t o of the diagonal coefficient Ru, where i is either 1 or 2, which is indicated by a 
horizontal dashed line on a figure. To find the value of an for which Ru, pc = Ru^o we draw a perpendicular from the 
point where it crosses the dotted line to the bottom axes. To find the value of <Ji for which Ru^in = Ru, w ^ draw 
a perpendicular from the point where the horizontal dashed line crosses the solid line to the top axes. For instance, 
the value i?22,o = 1-1 corresponds to an = 3 and a% = 0.62. The value an = 3, in its turn, gives i?n,o = 1-1 which 
corresponds to o\ — 0.54. 

One may start by specifying an, rather then Ru^q, to find a\ and 172. Then we draw a perpendicular from the 
bottom axes until it crosses the dotted line, which gives the value of Ru_ pc . Given the value of Ru^in to be the same, 
we find the value of Oi as described above. For the above example an = 3 corresponds to a\ = 0.54 and 02 = 0.62. 
We see, that we may not specify both o\ and &i independently: they must have the values which both correspond to 
the same an. For similar components, like those we are interested in, o\ and 02 should not differ much from each 
other, and therefore an, a coefficient which is related to the diffusion of one component through the other, should 
reflect this difference. 
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Having the diagonal coefficient mapped we have the parameters a qq and an defined uniquely (and taking into 
account that a\ q has negligible effect), as well as a\ and 172 fo r kinetic theory. We now compare the values of the 
cross coefficients given by "perturbation cell" method and kinetic theory. 



TABLE I: Gas-side transport coefficients obtained from kinetic theory 
and by "perturbation cell" method at T e? = 330 and fJ.i2, eq = 700 for 
P = 0.0002. 



parameters 


Rqq 


Ru 


R22 


Rql 


Rq2 


R12 


cri = 0.54 
cr 2 = 0.62 


2.96792e-011 


1. 11091 


1.09136 


3.82826e-007 


4.41483e-007 


0.0130511 


Otqq — 9 














Ctlq = 


3.01874e-011 


1.12461 


1.13991 


2.31477e-006 


2.27003e-006 


-0.816559 


an = 3 















One can see from Table [I] that while the diagonal coefficients are the sam^f], the cross coefficients we find are an 
order of magnitude larger than those found by kinetic theory. R12 even has a different sign. 

B. Temperature and chemical potential difference dependence 

In this subsection we investigate the dependence of the resistivity coefficients on the temperature and the chemical 
potential difference. On Fig. [SJS] we plot the these dependencies for R qq , R ql and Ru coefficients obtained from 
kinetic theory and "perturbation cell" method for the range of temperatures [325, . . . , 335] and for the range of 
chemical potential differences [400, . . . , 1000]. 




FIG. 3. Dependence of R qq on T and fii2 obtained from kinetic theory for ai = 0.54 and 02 = 0.62 (plane) and by "perturbation 
cell" method for a qq — 9, ai q = and an = 3 (points). 

The domain of T and /X12 is not big, so the dependence on them is linear, as expected. 



C. Validity of integral relations 

We compare the resistances found from numerical procedure to the values obtained from Eq. (IV.ll) . The relative 
difference between them is almost the same within the range of temperatures and chemical potential differences 



One should not expect exact compatibility between kinetic theory, which is most appropriate for gases with short range potentials, and 
the gradient theory, which is most appropriate for fluids with long range potentials. The purpose of this comparison in not to determine 
the exact values of adjustable parameters, but to show that it is possible to match coefficients in the two theories and to show the 
typical values of the parameters. 
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T/K 325 4 00 

1^2/J/mol 

FIG. 4. Dependence of R q i on T and /ii2 obtained from kinetic theory for o\ = 0.54 and 02 = 0.62 (plane) and by "perturbation 
cell" method for a qq — 9, ai q — and an = 3 (points). 



1.3> 




T/K 325 400 u^/J/mol 



FIG. 5. Dependence of Rn on T and /ii2 obtained from kinetic theory for o\ = 0.54 and 02 = 0.62 (plane) and by "perturbation 
cell" method for at qq = 9, a\ q = and Qn = 3 (points). 

considered: T — {325, • • • , 335} and /112 = {400, • • ■ , 1000}. In Table p] we give the relative errors for the resistances 
both for the case that we use the measurable heat fluxes on the vapor and on the liquid side. We refer to [35[ for 
details of the definition of the resistances using the measurable heat flux on the liquid side. 

TABLE II: Relative error in percent between the gas- and liquid- side co- 
efficients obtained by "perturbation cell" and "integral relations" meth- 
ods at T eq = 330 and fxi 2 , eq = 700 for /3 = 0.0002 and a qq = 9, ai q = 0, 
an = 3 . 

phase R qq Rn ^22 R^i Rg2 ^12 

gas 0.019090 0.064642 0.058851 0.020649 0.020680 0.097096 
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liquid 0.019090 0.006266 0.000432 0.036270 0.034886 6.233983 
The relative differences are not more then a few promille. It is larger only for R\ 2 which is discussed below. 

VII. DISCUSSION AND CONCLUSIONS 

In this paper we have studied stationary transport of heat and mass through the liquid- vapor interface in a mixture 



We used the expression for the excess entropy production of a surface derived from the continuous description 
which is identical to the one derived directly for the discrete description using the property of local equilibrium 
This makes it possible to give the linear force-flux relations for this case. These relations involve the interfacial 
resistances, which were the main focus of interest in this paper. Given the numerical solutions of the non-equilibrium 
square gradient model we were able to calculate these coefficients directly for a two-component mixture. Furthermore, 
we calculated these coefficients using integral relations, derived in [36| . This gives an independent way to determine 
the interfacial resistances. 

The main input parameters of the model are the local resistivity profiles used to calculate the continuous solution. 
There is not much theoretical information about the numerical value of these resistivities. In the vapor phase one can 
use kinetic theory. In the liquid phase it is most appropriate to use experimental values. There is no experimental 
information about the local resistivities in the interfacial region. As the local resistivities change in the surface from 
one bulk value to the other, it is natural to assume that they contain a contribution similar to the profile of the order 
parameter. There is also evidence from molecular dynamics simulations for one-component systems [38| that there 
is a peak in the local resistivity in the surface. As we are in the framework of the gradient theory, it is naturally to 
assume that such peaks are caused by a square gradient term, which is similar to the gradient contribution to the 
Hclmholtz energy density in the interfacial region, namely |Vp| 2 . The amplitudes of these peaks are not given by any 
theory and were used as parameters. We therefore get that the three local resistivities for a two-component mixture 
have the form given in Eq. (|V.3|) . Thus we get three adjustable amplitudes, a qq , a\ q and an, two of which were 
found to contribute significantly to the value of the transfer coefficients. 

In order to determine the typical values of the a's we need to compare our results with independently obtained 
resistivities. Unfortunately, not much experimental data are available for multi-component resistivities and, to the 
best of our knowledge, no data are available for our system. Furthermore, no molecular dynamic simulations of 
these properties are available for mixtures. The only available source of comparison is kinetic theory, which gives the 
expressions for the interfacial resistivities or transfer coefficients given in Eq. (IIII.2I) . We therefore compare our results 
to kinetic theory. Having three adjustable parameters in the gradient theory, a qq , a\ q and an, and two adjustable 
parameters in kinetic theory, the condensation parameters o\ and 02 , we are able to match three diagonal coefficients 
Rqq, R\i and i?22- We found that R qq does not really depend on a\ q and an. This makes it possible to fit a qq using 
R qq alone. For the values of the temperature and chemical potentials considered this gave a qq ~ 9. We furthermore 
found that the interfacial resistivities did not really depend on a\ q . We therefore took this amplitude equal to zero. 
In kinetic theory Rn and R22 depend on the condensation coefficients a± and <y%, respectively. Choosing an = 3 
gives values for the condensation coefficients of 0.54 and 0.62. As the components considered are very similar it is to 
be expected that these coefficients are close to each other. The values of a's obtained from the matching are such 
that the excess entropy production of the surface is positive, the second law is obeyed and the Onsager relations are 
valid. Having found the values of the a's from the diagonal transfer coefficients the values of the cross coefficients 
follow. 

We found that the values of the cross coefficients, obtained by our method are an order of magnitude larger than 
those found from kinetic theory. This confirms results from molecular dynamics simulations [2!| for a one-component 
system, where it was found that increasing the range of the attractive potential increased in particular the cross 
coefficients substantially above the values predicted by kinetic theory. This is an interesting result, indicating that 
kinetic theory underestimates the transfer coefficients for real fluids. This also indicates, that the effect of coup ling 
will be important in the interfacial region. Experiments also confirm the importance of the cross coefficients |23l . l24j . 

The effect of cross coefficients can be related to the measurable quantities, such as measurable heat of transfer 
q* = ~R q i/ R qq . This quantity can be associated both with gas and liquid phases in accordance to the corresponding 
heat fluxes. The difference q*' 9 — q*' = ~{R qi — R qi )/R qq — ~{h\ ~ h l ieq ) is equal to the difference of partial 

enthalpies between gas and liquid in equilibriunQ. This quantity is substantial, which implies that q*' 9 — q*' e is 
also substantial. This, in turn, makes the difference between the cross coefficients on the vapor and the liquid side 



9 We note that R 9 qq = R qq 
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substantial. This gives a theoretical ground for the importance of coupling in the interfacial region. Experiments 
[23l |24| confirm the size and importance of the heat of transfer on the vapor side. 

We did the comparison for one value of the temperature and chemical potential only. If one extends the analysis 
to a larger domain, one finds that the a's depend on the temperature and the chemical potential difference; we refer 
to [HI in this context. The results of kinetic theory [30l - l33j and molecular dynamics [28[ both support the existence 
of a peak in the diagonal local resistivities and therefore the use of finite values for a qq and an- 

Furthermore, it was found, that the data obtained directly from non-equilibrium numerical solution agree with 
the ones obtained using integral relations, as is expected. This gives an alternative and easier way to determine 
non-equilibrium properties of the interfacial region, needing only equilibrium information about the system. In fact, 
as we speak of linear non-equilibrium thermodynamics, this is the way it should be. The interfacial resistances are 
determined from equilibrium properties, just like Green-Kubo relations involve only equilibrium information in order 
to determine the transport coefficients. 

Appendix A: Determining an optimal perturbation rate 

The value of the resistance R(T e(? , ^12, eq ) does not depend on the perturbation, given the perturbation is small 
enough. However, the magnitude fj of the perturbation which may be considered sufficiently small, has to be deter- 
mined empirically. This would require considering perturbations where fi is beyond the appropriate range and will 
make the empirical resistances R(T eg , /ii2, eq ) to be dependent on 6. 

In order to determine the appropriate range of perturbations, that is when fi is small enough to consider them 
linear, and at the same time, large enough, to not interfere with the accuracy of the numerical solution, we check the 
obtained resistances for the thermodynamic consistency. We have the following constraints, which they must obey 
for each T and ^12: 

- i) the cross coefficients of each R matrix must satisfy Onsager relations; 

- ii) the second law consistency; 

- iii) coefficients obtained on the gas and the liquid side of the surface must be related; 

We will use the first condition to determine the range of fi, while the two remaining will be used for the verification 
of the results obtained in the paper. 

1. Onsager reciprocal relations 

As shown by Onsager [3^] , the cross coefficients must be the same. We therefore have R q i = Ri q and Rji = Rij . 

We calculate the coefficients at the values of equilibrium temperature and chemical potential difference T eq = 330 
K and /112. eq = 700 J/mol for different values of the adjustable amplitudes a qq , ai q , and an- 

In Tables [IIltiTV] we give the relative error in percent for the gas-side cross coefficients \(Rfj — R^)/ Rfj\-100% as 
a function of fi for a qq = 0, a\ q — 0, an ~ obtained by different methods. 

TABLE III: Relative error in percent for gas-side cross-coefficients ob- 
tained by "perturbation cell" method at T eq = 330 and /112, eq = 700 for 
different fi and for a qq — 0, a\ q — 0, an = 0. 



fi 


Rql 


R q 2 


R12 


2.0e-002 


8.963066 


35.863259 


34.908631 


2.0e-003 


0.273286 


0.369082 


19.683274 


2.0e-004 


0.011726 


0.007231 


1.909391 


2.0e-005 


0.066375 


0.071266 


2.336652 


2.0e-006 


4.963895 


8.128243 


5.843913 



TABLE IV: Relative error in percent for gas-side cross-coefficients ob- 
tained by "experiment like" method at T eq — 330 and |Ui2, eg = 700 for 
different fi and for a qq = 0, ai q = 0, an = 0. 



fi 


Rql 


R q 2 


R12 


2.0e-002 


1.275105 


0.828600 


754.982200 


2.0e-003 


0.038759 


0.363715 


38.708981 


2.0e-004 


0.131868 


0.238584 


6.247648 


2.0e-005 


1.301483 


2.056102 


20.984734 


2.0e-006 


13.282959 


20.788752 


632.124504 
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As one can see, — 0.02 is really an extreme perturbation and the difference is rather large. When we decrease j3 to 
2e-4 the differences become small. As we further decrease /3 to 2e-6 the inaccuracy of the numerical solution become 
comparable to the size of the perturbation. We conclude that the values for j3 to 2e-4 are closest to the converged 
values and use them as such. 

In Tables [VllVI] we give the same data for the higher continuous resistivities with rather substantial peak, when 
a qq = 10, ctiq — 10 and an — 10. As one can see, the Onsager relations are fulfilled there again best for f3 — 2e-4 

TABLE V: Relative error in percent for gas-side cross-coefficients ob- 
tained by "perturbation cell" method at T eq — 330 and /112, eq = 700 for 
different ft and for a qq — 10, ct\ q — 10, an = 10. 



p 


Rql 


R Q 2 


R12 


2.0e-002 


71.515410 


78.166809 


23.572836 


2.0e-003 


0.745604 


0.896547 


0.317348 


2.0e-004 


0.012358 


0.012650 


0.001919 


2.0e-005 


0.012078 


0.007485 


0.005290 


2.0e-006 


0.713969 


1.124994 


0.022121 



TABLE VI: Relative error in percent for gas-side cross-coefficients ob- 



tained by "experiment like" method at T eq = 330 and /U12, eq 
different /3 and for a qq = 10, ai q = 10, an = 10. 



700 for 



P 



Rql 



Rq2 



Rl 



2.0e-002 


4.225362 


2.559393 


12.259260 


2.0e-003 


0.443944 


0.256804 


1.091842 


2.0e-004 


0.068621 


0.019788 


0.093041 


2.0e-005 


0.269764 


0.407090 


0.008844 


2.0e-006 


2.717575 


4.149484 


2.025054 



We may notice that the behavior of the resistivities with respect to (3 is independent on the behavior of the 
resistivities with respect to a qq , ot\ q and an. This is natural, as these parameters control the different aspects of the 
system: /3 controls the perturbation rate, while a's are adjustable parameters, which control the size of the peak in 
the continuous resistivities. 



Appendix B: Consistency of the non-equilibrium solution 

1. Second law consistency 

In this subsection we investigate the values of parameters a qq , a,\ q , an for which the second law of thermodynamics 
is fulfilled. That is that the excess entropy production is positive and therefore the matrix of the resistivity coefficients 
is positive definite. This requires that the diagonal coefficients are positive and for each pair ql, q2 and 12 of the 
cross coefficients the expression 

DRik = RuRkk — -^{Rik + Rki) 2 > (B.l) 

must be positive. 

In Table |VII| we give the diagonal coefficients and expression (IB.4|) for each pair of the cross coefficients as a 
function of a qq for ai q — 0, an = and j3 — 2e-4 obtained by the "perturbation cell" method. In Tables [VIIHIIX] 
we give the same quantities for other choices of a. 

TABLE VII: 2nd law consistency for gas-side coefficients. The diagonal 
coefficients and the quantities defined by (|B.4|) . Data are obtained by 
"perturbation cell" method at T eq = 330 and /U12, eg = 700 for different 
a qq and for p = 0.0002, ai q = 0, an = 0. 



a qq 


Rqq 


Rn 


R22 


DR q i 


DR q2 


DR12 





7.05644e-015 


0.0754717 


-0.0919278 


2.13025e-015 


-2.59473e-015 


-0.0277518 


1 


3.36047e-012 


0.0937784 


-0.0741586 


1.26056e-012 


-9.9683e-013 


-0.0278179 


10 


3.35408e-011 


0.259425 


0.0851534 


3.48053e-011 


1.14244e-011 


0.0874467 
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TABLE VIII: 2nd law consistency for gas-side coefficients. The diagonal 
coefficients and the quantities defined by (|B.4[) . Data are obtained by 
"perturbation cell" method at T eq — 330 and /U12, eg = 700 for different 
ceiq and for (3 = 0.0002, a qq — 0, an = 0. 



Otlq 


Rqq 


i?n 


R22 


DR ql 


DR q2 


DR 12 





7.05644e-015 


0.0754717 


-0.0919278 


2.13025e-015 


-2.59473e-015 


-0.0277518 


1 


7.05608e-015 


0.0746391 


-0.0910331 


2.10664e-015 


-2.56935e-015 


-0.0271785 


10 


7.05304e-015 


0.0670813 


-0.0828915 


1.89251e-015 


-2.33855e-015 


-0.0222419 



TABLE IX: 2nd law consistency for gas-side coefficients. The diagonal 
coefficients and the quantities defined by (|B.4[) . Data are obtained by 
"perturbation cell" method at T eq — 330 and /X12, eg = 700 for different 
an and for /3 = 0.0002, a qq = 0, a lq = 0. 



an 


Rqq 


Rn 


R22 


DR ql 


DR q2 


DR 12 





7.05644e-015 


0.0754717 


-0.0919278 


2.13025e-015 


-2.59473e-015 


-0.0277518 


1 


7.05717e-015 


0.370078 


0.265626 


1.04468e-014 


7.49827e-015 


0.381226 


10 


7.10378e-015 


3.02063 


3.48284 


8.58316e-014 


9.89654e-014 


-69.2846 



We see, that the required quantities become positive for rather big values of a qq . They almost do not depend on 
the value of U\ q and they are positive for moderate values of parameter a\\. It is clear that finite values of a qq and 
an are needed to have a positive excess entropy production. 

All the above quantities almost do not depend on the value of j3 in the range [le-5, le-3]. The "experimental-like" 
procedure leads to almost the same values of all the quantities. The liquid-side coefficients reveal a similar behavior. 

2. Gas- and liquid- coefficients 

One can use the measurable flux J' q extrapolated from the liquid side of the surface, rather then from the gas side, 
using Eq. (|II.2|) . In this case one should use not the enthalpy of the gas bulk h 9 , but enthalpy of the liquid bulk h 
extrapolated to the diving surface. The two measurable fluxes are related as 

n 
i=l 

Identifying the forces and fluxes and writing the linear force-flux relations for the measurable heat flux on the liquid 
side, one introduces the interfacial resistances measured on the liquid side in the same way as it was done in Sec. [UJ 
for the interfacial resistances measured on the gas side. These resistances are related as follows 

R e = R 9 

Q1 11 

T>i _i_ 1 £ td£ _ r>5 _i_ u9 r>g 
^qi ' i, eq qq qi ' i, eq qq 

tdI 1 ui r>£ . j-? 9 -i— /i 9 jyg (2-3) 

ll iq > 'H, eq lx qq ~ ll iq > 'H, eq lx qq 

P! 1 4- h l R l 4- h l I?t 4- h& h& 7?^ - Ti 9 4- h 9 U 9 4- h 9 U 9 4- h 9 h 9 f/9 
"■ji ' 'H, eq ^jq ' n j, eq ^qi ' 'H, eq 11 j, eq ^qq ~ n ji ' r H, eq 1X1 jq ' 11 j, eq ^qi ' "i, eq IL j, eq ^qq 

These coefficients can be calculated independently from a non-equilibrium numerical solution. Given that, the validity 
of Eq. (|B.3p would indicate the internal consistency of the model. In this subsection we verify these relations. 
In Table [X] we give the relative error in percent between the left hand side and the right hand side of Eq. (|B.3() . 

TABLE X: Relative error in percent for invariant expressions in Eq. (|B.3|) 
obtained by "perturbation cell" method at T eq — 330 and Hi 2 , eq = 700 
for ft = 0.0002 and a qq — 1 , a lq = 1, an = 1. 

qq 11 22 ql lq q2 2q 12 21 
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0.000000 0.000002 0.000085 0.000001 0.000389 0.000001 0.000389 0.000060 0.000003 

For instance, the gl quantity is equal to \(R ql — h{ eq R qq ) — (R ql — h{ eq R qq )\/\R q i — h[ eq R e qq \-100%. The other 
quantities are defined in the same way. These errors almost do not depend neither on the value ol (3 in the range 
[le-5, le-3] nor on the values of a qq , ot\ q , a\\. The "experimental-like" procedure leads to almost the same results. 

3. Integral relations 

For two component mixture the force-flux equations have a form 

X q = R' qq J q — R ql — R q2 J% 2 

Xl = R{ q J q - Ru J^x - R \2 % (B.4) 
Xl = R'2q Jq — R2I J£.l — R22 

The left hand side of each equation must be equal to the right hand side. The difference therefore reflects the error. 
We give the relative error between the left an the right hand side of Eq. (|B.4[) in percent in Table |XIj . As a testing 
perturbation we used one of those used in the perturbation cell method. 

TABLE XI: Relative error in percent between the left- and right- hand 
side of Eq. (|B.4|) for coefficients obtained by "perturbation cell" and 
"integral relations" methods at T eq = 330 and fii2,e q = 700 for /3 — 
0.0002 and a qq = 9, ai q — 0, an = 3 . 







Integral relations 




Perturbation cell 


phase 


x q 


Xi X2 


x q 


Xi X2 


gas 


0.059489 


0.037918 0.296959 


0.046965 


0.087411 0.867098 


liquid 


0.059489 


0.172608 0.027275 


0.046851 


0.216819 0.014248 



Again, the relative difference is not more then a few promille. Given that this is the case even for a few percent 
difference in one of the coefficients, we may conclude that the values of the forces are insensitive to the precise value 
of this resistivity coefficient. This also indicates that the value of this coefficient obtained in [Z(| has a 6% error. This 
does not necessarily affect, however, the accuracy of the integral relations. 
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